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SCOPE  AND  PURPOSE 


It  is  typically  very  difficult  to  obtain  waiting-time  distribution 
functions  for  non-Markovian,  general  single-server  queues.  It  has  long  been 
recognized  that  virtually  any  probability  density  (such  as  for  service  times) 
can  be  approximated  quite  accurately  by  a  linear  combination  of  exponential 
density  functions.  Fortunately,  when  such  forms  can  be  assumed  to  describe 
interarrival  and  service  times,  classical  results  from  the  theory  of  queues 
permit  the  solution  of  the  waiting-time  problem  for  the  general  queue. 
However,  there  have  been  no  good  computer-based  algorithms  available  until 
recently  for  determining  appropriate  mixed  exponential  densities  to  fit  data 
or  closely  approximate  another  non-standard  density.  The  subject  of  this  work 
then  is  the  adaptation  of  a  special  nonlinear  optimization  routine  for 
generalized  mixture  estimation  and  its  use  in  calculating  very  accurate 
approximations  for  the  delay  distribution  of  any  general  single-server  queue. 
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I .  INTRODUCTION 
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In  a  recent  paper,  Fredericks  (1982)  offered  a  class  of  approximations 
for  the  GI/G/1  stationary  waiting-time  distribution  function  (CDF),  W^(t). 
The  main  idea  there  was  to  assume  a  specific  form  for  W^(t)  and  then  to  use 
the  Lindley  integral  equation  to  obtain  estimates  of  the  parameters. 

The  primary  form  used  for  the  approximant  was  the  exponential 

W  (t)  =  1  -  Ce"at 

q 

and  then  numerous  ways  were  documented  for  the  estimation  of  C  and  a.  How¬ 
ever,  presetting  the  form  of  W^(t)  leads  to  complications  when  more  complex 
functional  forms  are  tried,  and  this  approach  does  not  generally  portray  a 
waiting-time  CDF  in  an  accurate  way  over  its  full  range,  Fredericks  also 
mentioned  the  mixed  exponential  as  another  possible  form  for  W^(t)  in  light  of 
its  common  appearance  as  a  waiting-time  distribution.1  This  gave  a  more 
precise  approximation,  but  other  problems  then  surfaced. 

However,  there  is  an  entirely  different  way  to  look  at  the  problem  which 
has  some  precedence  in  the  literature,  and  which  often  provides  both  a  special 
insight  into  the  underlying  queue  mechanisms  and  a  most  accurate  way  of 
assessing  the  waiting  times.  We  suggest  looking  instead  at  the  interarrival 
and  service-time  distributions  defining  the  GI/G/1  as  the  candidates  for  the 
approximation,  such  that  the  final  form  of  W^(t)  is  a  member  of  a  very 


H 


s 
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By  mixed  exponential,  we  mean  that  the  complementary  CDF  1-W  (t)  may  be 


written  as  a  linear  combination  of  exponential  functions  with  possibly  complex 
powers.  In  its  most  general  form,  the  (possibly  complex)  linear  combination 
need  not  be  convex  and  the  constants  could  be  mixed  positive  and  negative. 
And,  of  course,  the  exponents  should  have  negative  real  parts.  See  Dehon  and 
Latouche  (1982)  for  an  expanded  and  up-to-date  discussion  of  this  class  of 
distribution  functions. 
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comprehensive  family  of  CDFs,  namely,  the  mixed  exponentials.  So  we  present  a 
way  to  approximate  which  guarantees  that  the  waits  turn  out  to  be  mixed 
exponentially  distributed. 


The  primary  motivation  for  our  alternative  approach  is  previous  work  of 
Smith  (1953),  and  Marchal  and  Harris  (1976).  The  major  relevant  result  of 
Smith  was  his  Theorem  4,  in  which  some  quite  general  sufficient  conditions 
were  presented  for  the  GI/G/1  waiting  times  to  have  a  mixed  exponential 
distribution.  Marchal  and  Harris  built  on  this  result  and  some  others  of 
Smith  to  offer  a  relatively  simple  mixed-exponential  approximation  of  W^(t) 
derived  by  fitting  the  difference  of  the  service  and  interarrival  time  random 
variables  by  a  difference  of  two  Erlang  variables. 


II.  NOTATION  AND  BACKGROUND  MATERIAL 


The  fundamental  results  for  the  GI/G/1  may  be  found,  for  example,  in 
Gross  and  Harris  (1974),  Chapter  6.  It  was  noted  there  that  the  GI/G/1 
problem  can  be  greatly  simplified  if  it  may  be  assumed  that  both  interarrival 
and  service  distributions  are  generalized  Erlangs  (GE)  expressed  as  con¬ 
volutions  of  independent  and  not-necessarily-identical  exponential  random 
variables.  When  the  means  of  such  exponentials  are  allowed  to  come  in  conju¬ 
gate  pairs  (so  that  their  Laplace-Stielt jes  transforms  are  inverse  poly¬ 
nomials),  Smith  calls  the  family  (with  n  the  degree  of  the  defining  poly¬ 
nomial).  Other  authors  (e.g.,  Cohen,  1982)  define  K  as  the  class  of  dis- 

n 

tributions  whose  transforms  are  rational  functions  (clearly  including  the 

inverse  polynomials);  but  we  shall  instead  call  these  R  (with  n  the  degree  of 

n 

the  denominator's  polynomial).  The  class  includes  all  regular  Erlangs,  but 

not  all  mixed  exponentials  and  mixed  Erlangs,  which  are  however  members  of  R^ . 

We  also  mention  the  generalized  phase-type  distributions  (PH)  popularized  by 

Neuts  and  others  (see  Neuts,  1981),  which  have  rational  transforms  as  well, 

though  not  necessarily  of  the  inverse  polynomial  form.  Thus  we  may 

symbolically  represent  the  relationship  of  those  respective  families  as 

GEC  K  C  R  and  PHC  R  • 
n  n  ^  n 

Now,  under  a  double  assumption  (i.e.,  that  the  queue  is  K^/K^/l),  it 
turns  out  that  the  delay  CDF  is  (for  example,  see  Gross  and  Harris) 

n  z  .t 

W  (t)  =  1  +  I  k.e  1  ,  (1) 

q  i=l  1 

where  the  {k^}  and  {z^}  would  be  determined  in  the  usual  way  (as  in  Gross  and 
Harris).  The  {k^}  are  arbitrary  in  sign,  while  the  {z^}  have  negative  real 
parts.  A  completely  parallel  result  exists  for  distributions  in  R  . 


Importantly,  Smith  also  showed  that  a  comparable  result  follows  even  for 
arbitrary  interarrival  times.  That  is,  the  GI/K^/1  and  GI/Rn/l  queues  have 
mixed  exponential  waits  independent  of  the  form  of  GI.  And  therein  lies  the 
key  to  our  approximation  method:  the  interarrival  and  service  CDFs  are  the 
functions  of  concern. 

So  our  approach  is  built  around  the  approximation  of  the  service  dis¬ 
tribution  by  a  member  of  either  the  class  K  or  R  .  Since  R  is  the  more 
J  n  n  n 

general  class  and  includes  K  ,  we  focus  on  selecting  from  amongst  its  members. 

n 

More  precisely,  we  work  with  a  subset  of  the  class  R^  made  up  of  those 

rational  functions  whose  denominators  have  real  roots,  and  call  this  class  GH 

2 

for  "generalized  hyperexponential." 

It  is  most  important  to  recognize  that  GH  is  a  very  complete  set  of 

potential  approximants  even  without  the  use  of  complex  scale  parameters.  The 

most  critical  characterization  of  this  coverage  is  the  fact  that  all  functions 

in  L2(0,«)  can  be  approximated  arbitrarily  closely  by  a  finite  linear  com- 

•At  +  3 

bination  of  functions  of  the  form  e  ,  $eR  .  (For  example,  see  Naylor  and 
Sell,  1971,  for  a  discussion  of  this  and  related  problems  on  the  Hilbert  space 
of  square  integrable  functions.) 

So  the  distribution  selection  problem  is  that  of  describing  the  proba¬ 
bility  structure  of  the  service  times  by  a  density  of  the  form 

K  -0 .  t 

b(t)  =  I  *.*.  e 
i=l 

2 

Note  that  the  class  GH  35GE  and  that  GHCPH. 

3 

Indeed,  the  classes,  GH,  PH  and  R  are  each  dense  in  the  set  of  all  CDFs 

n 

on  the  nonnegative  reals. 


either  from  data  or  as  an  approximation  to  an  actual  available  density 
function  whose  form  is  awkward.  If  we  indeed  have  raw  data,  then  our  problem 
is  just  one  of  parameter  estimation.  Otherwise,  we  have  a  curve  fitting  prob¬ 
lem,  which  in  turn  can  be  converted  to  an  estimation  problem  by  picking  a 
large  number  of  the  function's  points  and  then  proceeding  with  the  maximum- 
likelihood  estimation  (MLE).  The  recommended  approach  for  this  selection  is 
based  on  a  nonlinear  optimization  routine  previously  developed  for  the  MLE  of 
mixed  Weibull  parameters.  It  is  presented  in  the  next  section. 


III.  METHOD  FOR  SERVICE  DENSITY  APPROXIMATION 


The  numerical  procedure  for  estimation  has  been  built  up  from  previous 
work  on  exponential  and  Weibull  mixtures.  To  illustrate,  let  us  assume  that 
the  data  sampling  is  complete  so  that  all  service  times  are  fully  observed. 
In  the  event  that  there  are  incomplete  data  observations,  the  algorithm  is 
easily  altered. 

Maximum- likelihood  estimation  is  the  method  selected  mainly  because, 
under  fairly  general  conditions,  it  enjoys  the  important  limiting  statistical 
properties  of  efficiency,  normality,  and  unbiasedness.  Furthermore,  the  MLEs 
are  consistent,  invariant,  and  are  functions  of  sufficient  statistics  if  they 
exist.  When  sufficiency  and  unbiasedness  both  hold,  the  MLEs  are  also  of 
minimum  variance. 

A  first  key  observation  is  that  it  is  not  possible  to  obtain  explicit 
formulas  for  the  maximum- likelihood  estimators  of  parameters  involved  in  mixed 
exponential  densities  by  taking  the  partial  derivatives  and  equating  them  to 
zero.  Hence  we  must  resort  to  other  optimization  methods  and  numerical 
techniques.  Furthermore,  we  need  to  take  into  account  a  set  of  constraints  in 
addition  to  the  objective  function.  The  mixing  proportions  and  scale 
parameters  must  satisfy  some  simple  linear  relationships  and  there  may  exist 
other  constraints  related  to  the  sub-population  parameters.  Note  that  the 
constraints  are  generally  of  a  linear  type;  hence  the  problem  can  be  described 
as  a  mathematical  program  with  a  nonlinear  objective  function  and  linear 
constraints . 

The  target  criterion  function  in  the  maximum- likelihood  problem  is  the 
joint  density  function  for  a  random  sample  from  the  (mixed)  population 
governed  by  b(t).  As  is  <  'mmon,  it  is  much  easier  in  this  situation  to  work 


with  the  logarithm  of  the  likelihood  function.  Thus  if  we  write  the  likeli¬ 
hood  for  the  random  sample  t^, . . . ,t^  as 

N  N  K  -<t>.  t. 

L(a)  =  n  b(t.;a)  =  II  l  *.*.e  3  1  (5) 

i=l  1  i=l  j=l  J  J 

where  a  is  the  vector  of  parameters  (which  may  include  K) ,  then  its  logarithm 
is  expressed  as 

N 

£(a)  =  £  In  b(t.;a).  (6) 

i=l  1 

The  general  MLE  problem  for  the  mixture  may  then  be  formulated  as  the  non¬ 
linear  constrained  optimization  problem: 

max  £(a) 

a 

subject  to 

a  t  S  =  {a|£2L  =  1;*  >  0).  (7) 

Under  the  standard  mixed-exponential  regime,  all  fr.  >  0,  and  .  would  be 

J  J 

real  and  greater  than  0.  The  most  efficient  algorithm  available  for  the 
solution  of  this  problem  is  due  to  the  joint  efforts  of  Kaylan  (1978)  and 
Kaylan  and  Harris  (1981),  and  Mandelbaum  (1982)  and  Mandelbaum  and  Harris 

(1982).  It  is  an  iterative  numerical  procedure  built  around  Newton's  method 

with  additional  use  of  second  partial  derivatives  to  speed  up  convergence 
whenever  necessary. 

Because  the  mixing  parameters  may  be  negative  in  our  queueing  problem, 
the  Kaylan/Mandelbaum/Harris  algorithm  has  been  carefully  altered.  The  basic 
approach  is  similar,  except  that  now  the  number  of  mixing  variables  has  been 
doubled.  Since  it  is  preferable  to  optimize  over  nonnegative  variables,  we 
thus  write  JT  as  the  difference  between  its  positive  and  negative  parts. 

There  were  two  major  changes  necessary  in  the  algorithm.  First,  additional 


code  had  to  be  added  to  make  sure  that  the  density  function  b(.t)  did  not 
become  negative.  A  second  alteration  had  to  be  made  to  prevent  the  mixing 
parameters  {S\}  from  drifting  too  far  out  along  the  real  line  in  either 
direction. 


IV.  ILLUSTRATIVE  EXAMPLE 


To  show  how  our  proposed  method  would  work,  the  following  test  problem  is 
offered.  It  is  patterned  somewhat  after  an  example  in  Gross  and  Harris 
(pp.  300-301)  and  is  of  the  M/G/l  type.  Interarrivals  are  assumed  to  have 
mean  2,  while  a  set  of  50  service  times  has  been  generated  randomly  according 
to  an  Erlang(2)  with  mean  1.  These  numbers  are  (to  two  decimal  places): 


.90, 

2.14, 

1.33, 

.52, 

•  72, 

.06, 

•  49, 

■53, 

1.10, 

.63 

•31, 

1.47, 

.91, 

.26, 

.49, 

.97, 

.24, 

2.99, 

2.60, 

1.11 

1.99, 

.41, 

31, 

1.58, 

1.42, 

.84, 

.11, 

•  54. 

1.27, 

2.05 

.51, 

.91, 

1.78, 

1.69, 

.78, 

.64, 

1.62, 

1.96, 

.93, 

.41 

1.70, 

.73, 

.98, 

.71, 

.73, 

.44, 

1.27, 

.92, 

.71, 

.78 

The  algorithm  of  the  previous  section  was  then  applied  to  the  data  and 
the  resultant  generalized  mixed  exponential  turned  out  to  be  (with  its 
parameters  rounded  off  for  simplicity) 

B(t)  =  1  -  e  +  4e  -  4e 

The  next  step  is  to  substitute  B(t)  into  the  GI/G/1  waiting-time 
equations  as  those  presented  earlier  in  Gross  and  Harris.  It  follows  that 
W^(t)  may  be  written  as 

W  (C)  =  l  +  t  .-8779t  +  k  -<3.811+5335i)t  +  -<3.8U-.5335i)t 

q  1  2  3 

In  order  to  guarantee  real  values  for  W^(t),  k^  must  be  the  complex 
conjugate  of  k^.  The  result  can  then  be  simplified  to 

Vq  ( t )  =  1  +  k^e"' 8779t  1  k2,-(J.811+.53351)t  +  ^-(3.611- .  5335i)t 

When  the  work  of  Smith  is  applied  and  the  corresponding  residues  obtained,  we 
find  that  k^  =  -.5141  and  =  .02884  +  .03575i.  The  final  expression  for  the 
stationary  waiting-time  CDF  is  thus 

Wq(t)  =  1  -  ,5141e',8779t  +  e‘3-811t[ .05768  cos  ,5335t  -  .07150  sin  .5335t]. 


V.  CONCLUDING  REMARKS 


It  should  be  noted  that  several  other  methods  are  available  in  the 
literature  which  provide  similar  approximating  forms  for  the  service-time 
distribution.  In  particular,  Luchak  (1956)  and  Wishart  (1959)  expanded  the 
Stieltjes  transform  B*(s)  in  a  Laurent  series  about  (-y),  and  then  truncated 
after  k  terms  to  give 

k 

B*(s)  =  Z  c  ( 1+s/y) ~m  . 

_ .  m 
m=l 

(This  corresponds  to  approximating  B(t)  as  a  mixture  of  Erlangs  with  the  same 

scale.)  But  these  {c  }  are  not  likely  to  sum  to  one,  and  thus  the  truncation 

m 

is  not  a  true  CDF.  We  may  try  instead  to  match  moments  and  use  boundary 
conditions  for  estimating  y  and  the  {cm}.  Unfortunately,  the  latter  estimates 
may  not  be  the  same  as  those  from  the  Laurent  expansion  and  would  not  neces¬ 
sarily  possess  any  of  the  major  properties  typically  desired.  Observe  that 

such  Laurent  forms  are  indeed  members  of  the  class  R  . 

n 

In  closing  now,  we  believe  it  important  to  view  this  work  as  a  total 
package:  numerical,  statistical  and  mathematical.  We  see  our  approach  as  a 
unification  of  some  of  the  most  fundamental  results  on  the  GI/G/1  queue, 
particularly  building  on  the  early  and  basic  work  of  Smith.  This  and  related 
work  which  have  followed  are  all  put  in  total  perspective  by  Cohen.  The 
numerical  portions  of  this  study  began  with  the  optimization  procedure  and 
ultimately  led  to  the  derivation  of  the  stationary  delay  CDF.  But  we  really 
cannot  separate  this  from  the  attempt  to  find  quality  statistical  estimators. 
By  virtue  of  the  fact  that  the  distribution  selection  is  performed  by  MLE,  it 
is  clear  that  this  sort  of  complete  modeling  effort  is  very  much  data  oriented 
and  therefore  a  potentially  powerful  tool  for  the  applied  queueing  analyst. 


10 


Finally,  the  extreme  simplicity  of  the  mathematical  solution  to  the  waiting 
time  problem  should  be  reemphasized. 
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